home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DSINV.FOR < prev    next >
Text File  |  1985-11-29  |  4KB  |  123 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DSINV
  5. C
  6. C        PURPOSE
  7. C           INVERT A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
  8. C
  9. C        USAGE
  10. C           CALL DSINV(A,N,EPS,IER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           A      - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN
  14. C                    SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT
  15. C                    MATRIX.
  16. C                    ON RETURN A CONTAINS THE RESULTANT UPPER
  17. C                    TRIANGULAR MATRIX IN DOUBLE PRECISION.
  18. C           N      - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.
  19. C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED
  20. C                    AS RELATIVE TOLERANCE FOR TEST ON LOSS OF
  21. C                    SIGNIFICANCE.
  22. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  23. C                    IER=0  - NO ERROR
  24. C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
  25. C                             TER N OR BECAUSE SOME RADICAND IS NON-
  26. C                             POSITIVE (MATRIX A IS NOT POSITIVE
  27. C                             DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI-
  28. C                             FICANCE)
  29. C                    IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFI-
  30. C                             CANCE. THE RADICAND FORMED AT FACTORIZA-
  31. C                             TION STEP K+1 WAS STILL POSITIVE BUT NO
  32. C                             LONGER GREATER THAN ABS(EPS*A(K+1,K+1)).
  33. C
  34. C        REMARKS
  35. C           THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE
  36. C           STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS.
  37. C           IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU-
  38. C           LAR MATRIX IS STORED COLUMNWISE TOO.
  39. C           THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL
  40. C           CALCULATED RADICANDS ARE POSITIVE.
  41. C
  42. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  43. C           DMFSD
  44. C
  45. C        METHOD
  46. C           SOLUTION IS DONE USING FACTORIZATION BY SUBROUTINE DMFSD.
  47. C
  48. C     ..................................................................
  49. C
  50.       SUBROUTINE DSINV(A,N,EPS,IER)
  51. C
  52. C
  53.       DIMENSION A(1)
  54.       DOUBLE PRECISION A,DIN,WORK
  55. C
  56. C        FACTORIZE GIVEN MATRIX BY MEANS OF SUBROUTINE DMFSD
  57. C        A = TRANSPOSE(T) * T
  58.       CALL DMFSD(A,N,EPS,IER)
  59.       IF(IER) 9,1,1
  60. C
  61. C        INVERT UPPER TRIANGULAR MATRIX T
  62. C        PREPARE INVERSION-LOOP
  63.     1 IPIV=N*(N+1)/2
  64.       IND=IPIV
  65. C
  66. C        INITIALIZE INVERSION-LOOP
  67.       DO 6 I=1,N
  68.       DIN=1.D0/A(IPIV)
  69.       A(IPIV)=DIN
  70.       MIN=N
  71.       KEND=I-1
  72.       LANF=N-KEND
  73.       IF(KEND) 5,5,2
  74.     2 J=IND
  75. C
  76. C        INITIALIZE ROW-LOOP
  77.       DO 4 K=1,KEND
  78.       WORK=0.D0
  79.       MIN=MIN-1
  80.       LHOR=IPIV
  81.       LVER=J
  82. C
  83. C        START INNER LOOP
  84.       DO 3 L=LANF,MIN
  85.       LVER=LVER+1
  86.       LHOR=LHOR+L
  87.     3 WORK=WORK+A(LVER)*A(LHOR)
  88. C        END OF INNER LOOP
  89. C
  90.       A(J)=-WORK*DIN
  91.     4 J=J-MIN
  92. C        END OF ROW-LOOP
  93. C
  94.     5 IPIV=IPIV-MIN
  95.     6 IND=IND-1
  96. C        END OF INVERSION-LOOP
  97. C
  98. C        CALCULATE INVERSE(A) BY MEANS OF INVERSE(T)
  99. C        INVERSE(A) = INVERSE(T) * TRANSPOSE(INVERSE(T))
  100. C        INITIALIZE MULTIPLICATION-LOOP
  101.       DO 8 I=1,N
  102.       IPIV=IPIV+I
  103.       J=IPIV
  104. C
  105. C        INITIALIZE ROW-LOOP
  106.       DO 8 K=I,N
  107.       WORK=0.D0
  108.       LHOR=J
  109. C
  110. C        START INNER LOOP
  111.       DO 7 L=K,N
  112.       LVER=LHOR+K-I
  113.       WORK=WORK+A(LHOR)*A(LVER)
  114.     7 LHOR=LHOR+L
  115. C        END OF INNER LOOP
  116. C
  117.       A(J)=WORK
  118.     8 J=J+K
  119. C        END OF ROW- AND MULTIPLICATION-LOOP
  120. C
  121.     9 RETURN
  122.       END
  123.